Clasificacion_Supervisada
Jessica Bernal
20/3/2022
Técnicas de clasificación y evaluación de procesos en sistemas forestales
Jessica Bernal Borrego
Contenido:
- Introducción
- Creación muestreo
- Clasficación supervisada
- Clasificador KNN
- Clasificador ANN
- Clasificador SVM
- Clasificador RF
- comparativa resultados entre clasificadores
- Clasificador multi-temporal
###Introducción
Los objetivos de este tutorial son: - Clasificar una extracto de una escena Sentinel 2 usando diferentes clasificadores de machine learning. - Comparar los resultados obtenidos empleando distintos clasificadores
Los paquete que emplearemos son:
library(sp)
library(rgdal)
library(raster)
library(reshape)
library(grid)
library(gridExtra)
library(RStoolbox)
library(caret)
library(rasterVis)
library(corrplot)
library(doParallel)
library(NeuralNetTools)
library(tidyr)
library(stringr)
library(e1071)
library(sf)
library (mapview) El siguiente paso consistirá en definir el directorio de trabajo donde se localizarán nuestras imágenes. Los datos se corresponden con una escena Sentinel 2 (fichero denominado sentinel_o_bxx.tif, siendo xx el número de la banda)
setwd("C:/Geoforest/Tec_Clasif")
dir_in<-"./Material_practicas/Sentinel/O"Como en cuadernos anteriores creamos una lista con los nombres de los archivos alojados en el directorio de trabajo, creando posteriormente un rasterstack
rasList <- list.files(dir_in,pattern="tif",
full.names = TRUE)
sentinel_o <- stack(rasList)Una vez generado, vamos a proceder a comprobar los atributos de la escena
sentinel_o## class : RasterStack
## dimensions : 2052, 2057, 4220964, 10 (nrow, ncol, ncell, nlayers)
## resolution : 20, 20 (x, y)
## extent : 300920, 342060, 4042100, 4083140 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : sentinel_o_b01, sentinel_o_b02, sentinel_o_b03, sentinel_o_b04, sentinel_o_b05, sentinel_o_b06, sentinel_o_b07, sentinel_o_b08, sentinel_o_b09, sentinel_o_b10
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535A continuación se van a generar el gráfico de densidad y el histograma para una de las bandas. En primer lugar lo calcularemos para una de las bandas (imagen).
gdensidad=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_density()
ghistograma=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_histogram()
print (gdensidad)
print (ghistograma)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Si las queremos hacer de todas las bandas es posible generar un gráfico por cada una de ellas y posteriormente componerlos en una sola figura
gdens1=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_density()
gdens2=ggplot(sentinel_o,aes(sentinel_o_b02))+geom_density()
gdens3=ggplot(sentinel_o,aes(sentinel_o_b03))+geom_density()
gdens4=ggplot(sentinel_o,aes(sentinel_o_b04))+geom_density()
gdens5=ggplot(sentinel_o,aes(sentinel_o_b05))+geom_density()
gdens6=ggplot(sentinel_o,aes(sentinel_o_b06))+geom_density()
gdens7=ggplot(sentinel_o,aes(sentinel_o_b07))+geom_density()
gdens8=ggplot(sentinel_o,aes(sentinel_o_b08))+geom_density()
gdens9=ggplot(sentinel_o,aes(sentinel_o_b09))+geom_density()
gdens10=ggplot(sentinel_o,aes(sentinel_o_b10))+geom_density()
grid.arrange(gdens1,gdens2,gdens3,gdens4,gdens5,gdens6,gdens7,gdens8,gdens9,gdens10,ncol=4,nrow=4)
Creación muestreo
Al ser una clasificación supervisada necesitaremos aportar al clasificador la información necesaria para realizar las fases de entrenamiento y validación. A partir del MFE facilitado, con geometría poligonal, tenemos dos opciones a la hora de continuar. Para ello, es necesario saber que en este entorno de trabajo la información suministrada al clasificador deberá ser de tipo poligonal. Por ello, las opciones son:
Opción 1: Preparación de los datos desde un software externo, por ejemplo QGIS, y luego leerlo en R.
Opción 2: Preparación de los datos desde R.
En el caso de optar por trabajar con una herramienta externa los datos de tipo puntual pueden ser almacendados en formato shapefile y luego ser leidos mediante la función readOGR.
A modo de ejemplo, la llamada a la función sería:
#train_data<-readOGR('./entrenamiento/muestreo_500.shp')Debiendo comprobar posteriormente la estructura del fichero leido.
#str(train_data)Tambien es posible hacer la misma operación mediante la función st_read.
En este cuaderno se optará por trabajar con la segunda opción, es decir, preparando el muestreo desde R.
Para ello, en primer lugar se procederá a leer el archivo shapefile del MFE.
#Generamos una semilla para garantizar la repetitividad de los resultados
set.seed(123)
MFE=st_read('./Material_practicas/MFE/MFE.shp')## Reading layer `MFE' from data source
## `C:\Geoforest\Tec_Clasif\Material_practicas\MFE\MFE.shp' using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 310927.3 ymin: 4052098 xmax: 332056.3 ymax: 4073143
## Projected CRS: ETRS89 / UTM zone 30Nmapview (MFE,zcol='leyenda')Con objeto de obtener una muestra balanceada se va a determinar el total de la superficie ocupada por cada clase de la leyenda, repartiendo el tamaño de la muestra proporcional a la superficie ocupada.
clases = unique(MFE$leyenda)
area_total = sum(st_area(MFE))
area_clases=0
for (i in 1:length(clases)) {
geom_clase=MFE[which(MFE$leyenda == clases[i],arr.ind=FALSE),]
area_clases[i]=sum(st_area(geom_clase))
}Se va a fijar un tamaño total de muestreo igual a 500 puntos de tal forma que en cada clase habrá el siguiente número de muestras.
num_muestras= as.integer(500*area_clases/area_total)
print(num_muestras)## [1] 28 0 26 6 0 0 15 0 25 36 122 34 203En este punto, analizar el número de cada clase. ¿Hay muestras en todas las clases?, ¿Que significa que haya clases con un número elevado y otras que sea muy reducido o directamente cero?
Aun sabiendo que no es correcto, en lugar de establecer el muestreo atendiendo al criterio anterior se va a seleccionar un número fijo de muestras para todas las clases. De esta manera, el trabajo a entregar por el alumno será determinar una leyenda adecuada a la variabilidad espacial y espectral de las clases presentes en la escena.
Por ello, en primer lugar se va a proceder a realizar un muestro sobre el MFE de tipo aleatorio, extrayendo la información temática a partir de la función st_join.
puntos.ref <- st_sample(MFE, c(50,50,50,50,50,50,50,50,50,50,50,50,50), type='random',exact=TRUE) #Generamos una lista de puntos de forma aleatoria
puntos.ref<-st_sf(puntos.ref) #Convertimos la lista en un spatial feature
puntos.ref<-st_join(puntos.ref,MFE) #Cruzamos los datos
puntos.ref_backup <- puntos.ref
mapview(puntos.ref, zcol='leyenda')Y su representación sobre el MFE:
mapview(MFE, zcol='leyenda')+mapview(puntos.ref,zcol='leyenda')A continuación se va a proceder a obtener la firma espectral de cada una de las clases. en primer lugar será necesario extraer los valores de reflectancia para cada punto en cada una de las bandas mediante el comando extract. Como estos datos los representaremos mediante la librería ggplot los convertiremos a un tipo de dato dataframe. Además, antes de la representación se determinarán los valores medios de reflectancia por clase y banda.
Así, en primer lugar generamos el dataframe.
puntos.ref=as_Spatial(puntos.ref)
puntos.ref@data$leyenda=as.factor(puntos.ref@data$leyenda)
reflectancia<- as.data.frame(raster::extract(sentinel_o,puntos.ref))Comprobando los valores extraidos.
head(reflectancia)## sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05
## 1 935 701 434 684 1405
## 2 1069 871 703 985 2058
## 3 869 562 356 456 800
## 4 875 583 360 454 843
## 5 922 681 469 694 1312
## 6 929 696 471 702 1611
## sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10
## 1 1705 1545 1754 689 340
## 2 2270 2263 2531 1983 1174
## 3 962 787 883 497 258
## 4 853 901 1032 755 333
## 5 1720 1605 1697 1248 698
## 6 1998 1923 2179 1293 633Calculamos el valor medio de reflectancia de cada clase para cada banda. Para ello, usaremos la función aggregate() para unir los puntos de entrenamiento por clase.
mean_reflectancia <-aggregate(reflectancia,list(puntos.ref$leyenda),mean,na.rm = TRUE)Comprobamos los valores medios obtenidos
head(mean_reflectancia)## Group.1 sentinel_o_b01 sentinel_o_b02 sentinel_o_b03
## 1 -- 1088.940 954.56 782.28
## 2 Bosques Mixtos 1040.033 842.18 704.48
## 3 Cadufolios 1348.460 1171.52 1135.10
## 4 Castaños/Caducifolios 1105.320 954.24 1059.78
## 5 Eucalipto 926.980 672.34 475.60
## 6 Matorral 1460.080 1292.59 1219.90
## sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## 1 1126.10 2144.880 2429.600 2449.42 2619.32
## 2 910.86 1535.053 1722.407 1719.52 1846.22
## 3 1327.40 1725.660 1868.640 1904.24 2057.20
## 4 1233.26 1441.640 1580.820 1616.80 1753.70
## 5 629.52 1082.500 1208.240 1133.80 1279.34
## 6 1449.83 2027.840 2184.040 2173.93 2315.40
## sentinel_o_b09 sentinel_o_b10
## 1 1576.540 945.54
## 2 1255.907 784.16
## 3 1897.600 1326.12
## 4 2244.360 1609.96
## 5 784.000 439.12
## 6 1613.840 1137.70Por la forma en la que estan almacenados los datos en el dataframe (cada banda se almacena en una columna) es necesario modificarlo para que esten todos los datos de reflectancias registrados en una columna, creando una nueva columna donde se registre la banda de donde proceden, de tal forma que la información aparecerá ordenada por filas.
mean_reflectance2 <- gather(mean_reflectancia, key="banda", value="reflectance",sentinel_o_b01:sentinel_o_b10)Si analizamos el contenido del dataframe no se dispone de un campo numerico que permita ordenar las bandas a la hora de pintarlas. Por ello se va a crear un nuevo campo de tipo numérico con el número de la banda.
mean_reflectance2$banda_num=(as.numeric(str_replace(mean_reflectance2$banda,"sentinel_o_b","")))Finalmente, mediante ggplot se pintará un gráfico de tipo geom_line para ver la firma espectral de cada una de clases.
Como se puede ver, muchas de las clases presentan un comportamiento muy similar y por tanto la calidad temática de los resultados de la clasificación a priori pueden ser bajos.
ggplot(mean_reflectance2,aes(x=banda_num,y=reflectance))+
geom_line(aes(colour = Group.1))+theme_bw()
Clasficación supervisada
En primer lugar veremos el uso de alguno de los operadores clásicos empleados en Teledetección para clasificar imágenes, en este caso clasificador por máxima probabilidad. Para ello se empleará la función superClass dentro del paquete RStoolbox. De entre los parámetros a incluir en la función destacar que:
trainData contiene el muestreo a emplear en la clasificación.
trainPartition: contiene un valor indicando el tamaño de la muestra destinada al entrenamiento.
model indica el tipo de clasficación que se desea realizar, por defecto la función aplica un random forest (rf) pero en este caso se realizará una clasificación por máxima probabilidad (mlc).
#puntos.ref<- as_Spatial(puntos.ref)
puntos.ref@data$leyenda=as.factor(puntos.ref@data$leyenda)
Max.prob<- superClass(sentinel_o,
trainData = puntos.ref,
trainPartition =0.5,
responseCol = "leyenda",
model = "mlc") #máxima probabilidadEl resultado de la clasificación se muestra recogido en una variable de tipo lista. En el quinto elemento de la lista se recoge el resultado cartográfico mediante un rasterLayer, pudiendo ser representado por ejemplo mediante la función plot.
leyenda_colores <- viridis::viridis(13)
plot(Max.prob$map,
col=leyenda_colores,
legend = FALSE)
legend("topright",cex=0.65, y.intersp = 0.55,x.intersp = 0.5,
legend = levels(as.factor(puntos.ref$leyenda)),
fill = leyenda_colores ,title = "",
inset=c(0,0))
Además del producto cartográfico es necesario realizar una evaluación de la calidad temática. Para ello mediante en el segundo elemento de la lista, denominado modelFit se encuentra la matriz de confusión resultante así como los valores de exactitud global y kapp obtenidos en el entrenamiento. Por otro lado en el elemento results aparecen recogidos estos elementos de calidad global y su desviación.
Max.prob$modelFit## [[1]]
## TrainAccuracy TrainKappa method
## 1 0.354396 0.2797632 custom
##
## [[2]]
## Cross-Validated (5 fold) Confusion Matrix
##
## (entries are average cell counts across resamples)
##
## Reference
## Prediction -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 1.8 0.2 0.0 0.0
## Bosques Mixtos 0.2 3.8 1.0 0.0
## Cadufolios 0.0 0.0 1.0 0.4
## Castaños/Caducifolios 0.0 0.0 0.0 3.8
## Eucalipto 0.2 1.4 0.2 0.0
## Matorral 0.2 0.2 0.4 0.0
## Pinos 1.6 3.0 0.0 0.0
## Pinos/Pinsapos 0.2 0.8 0.0 0.0
## Pinsapos 0.0 0.4 0.0 0.0
## Quercineas 0.0 2.6 0.4 0.0
## Ribera 0.0 0.6 0.0 0.0
## Suelos? 0.2 0.8 1.0 0.2
## Veg. Dispersa 0.2 0.8 1.0 0.4
## Reference
## Prediction Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 0.0 0.0 1.2 0.6 0.0
## Bosques Mixtos 0.2 1.6 0.6 0.6 0.2
## Cadufolios 0.0 0.4 0.0 0.0 0.2
## Castaños/Caducifolios 0.0 0.0 0.0 0.0 0.0
## Eucalipto 1.8 0.6 0.6 0.4 0.2
## Matorral 0.2 3.6 0.0 0.0 0.2
## Pinos 0.6 0.8 10.8 1.0 1.0
## Pinos/Pinsapos 1.2 0.8 2.0 1.0 0.8
## Pinsapos 0.2 0.4 2.2 0.4 2.0
## Quercineas 0.6 0.4 1.4 1.0 0.0
## Ribera 0.0 0.0 0.2 0.0 0.2
## Suelos? 0.0 0.0 0.6 0.0 0.2
## Veg. Dispersa 0.2 1.4 0.4 0.0 0.0
## Reference
## Prediction Quercineas Ribera Suelos? Veg. Dispersa
## -- 0.2 0.0 0.2 0.2
## Bosques Mixtos 0.8 0.2 1.4 1.8
## Cadufolios 0.2 0.0 0.2 1.6
## Castaños/Caducifolios 0.0 0.0 0.0 0.0
## Eucalipto 0.8 0.2 0.0 0.4
## Matorral 0.0 0.2 0.0 1.6
## Pinos 0.8 1.2 1.4 1.6
## Pinos/Pinsapos 0.2 0.2 0.0 0.6
## Pinsapos 0.0 0.6 0.4 0.2
## Quercineas 6.6 0.8 0.0 4.2
## Ribera 0.4 0.2 0.0 0.0
## Suelos? 0.0 0.0 0.4 1.4
## Veg. Dispersa 0.0 0.2 1.0 1.4
##
## Accuracy (average) : 0.3544Max.prob$model$results## parameter Accuracy Kappa AccuracySD KappaSD
## 1 none 0.354396 0.2797632 0.02372322 0.02503495Creación de un dataframe con los puntos del entrenamiento etiquetados y con sus valores de reflectancia
En el paso anterior se generó un dataframe con los valores medios de reflectancia para cada clase. Ahora se va a generar un dataframe que contendrá para cada punto su etiqueta y los valores de reflectancia de todas y cada una de las bandas. Advertir que el objeto train_data@data apunta a toda esa información, de forma que *@* es un operador especial que permite acceder a un objeto dentro de otro objeto. Nota: Se va a crear una variable denominada train_data igual a puntos.ref por si el trascurso de la práctica se comete un error, pudiendo tener así un backup de los datos hasta este punto.
train_data = as_Spatial(puntos.ref_backup)
train_data@data$leyenda=as.factor(puntos.ref@data$leyenda)
train_data@data=data.frame(train_data@data,reflectancia[match(rownames(train_data@data),rownames(reflectancia)),])Veamos el resultado obtenido.
str(train_data)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 1100 obs. of 15 variables:
## .. ..$ OBJECTID : num [1:1100] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ NOM_FORARB : chr [1:1100] "Alcornocales" "Alcornocales" "Alcornocales" "Alcornocales" ...
## .. ..$ Shape_Leng : num [1:1100] 34046 34046 34046 34046 34046 ...
## .. ..$ Shape_Area : num [1:1100] 4419619 4419619 4419619 4419619 4419619 ...
## .. ..$ leyenda : Factor w/ 13 levels "--","Bosques Mixtos",..: 10 10 10 10 10 10 10 10 10 10 ...
## .. ..$ sentinel_o_b01: num [1:1100] 935 1069 869 875 922 ...
## .. ..$ sentinel_o_b02: num [1:1100] 701 871 562 583 681 696 609 753 500 654 ...
## .. ..$ sentinel_o_b03: num [1:1100] 434 703 356 360 469 471 377 542 282 391 ...
## .. ..$ sentinel_o_b04: num [1:1100] 684 985 456 454 694 702 477 791 254 562 ...
## .. ..$ sentinel_o_b05: num [1:1100] 1405 2058 800 843 1312 ...
## .. ..$ sentinel_o_b06: num [1:1100] 1705 2270 962 853 1720 ...
## .. ..$ sentinel_o_b07: num [1:1100] 1545 2263 787 901 1605 ...
## .. ..$ sentinel_o_b08: num [1:1100] 1754 2531 883 1032 1697 ...
## .. ..$ sentinel_o_b09: num [1:1100] 689 1983 497 755 1248 ...
## .. ..$ sentinel_o_b10: num [1:1100] 340 1174 258 333 698 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:1100, 1:2] 326741 329242 328773 331246 329865 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "lon" "lat"
## ..@ bbox : num [1:2, 1:2] 312560 4052227 331754 4072825
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "lon" "lat"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
## .. .. ..$ comment: chr "PROJCRS[\"ETRS89 / UTM zone 30N\",\n BASEGEOGCRS[\"ETRS89\",\n DATUM[\"European Terrestrial Reference"| __truncated__Podemos observar como nos indica que tenemos 11 variables, una correspondiente a la clase y diez a las bandas espectrales, siendo estas últimas nuestras variables predictoras de la variable respuesta consistente en las clases. Ahora, podriamos ver un resumen estadístico de la variable.
summary(train_data@data)## OBJECTID NOM_FORARB Shape_Leng Shape_Area
## Min. : 1.0 Length:1100 Min. : 1805 Min. : 52512
## 1st Qu.: 6.0 Class :character 1st Qu.: 3900 1st Qu.: 84213
## Median :11.5 Mode :character Median : 26445 Median : 2351264
## Mean :11.5 Mean : 64448 Mean : 9131237
## 3rd Qu.:17.0 3rd Qu.: 80456 3rd Qu.:10093586
## Max. :22.0 Max. :455519 Max. :81665992
##
## leyenda sentinel_o_b01 sentinel_o_b02 sentinel_o_b03
## Pinos :200 Min. : 761 Min. : 457.0 Min. : 249.0
## Bosques Mixtos:150 1st Qu.: 888 1st Qu.: 632.8 1st Qu.: 408.0
## Veg. Dispersa :150 Median : 980 Median : 787.5 Median : 587.5
## Matorral :100 Mean :1103 Mean : 905.1 Mean : 775.2
## Quercineas :100 3rd Qu.:1126 3rd Qu.: 995.2 3rd Qu.: 963.2
## -- : 50 Max. :5384 Max. :5290.0 Max. :5925.0
## (Other) :350
## sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## Min. : 226.0 Min. : 247 Min. : 233 Min. : 212 Min. : 187
## 1st Qu.: 561.5 1st Qu.:1030 1st Qu.:1158 1st Qu.:1146 1st Qu.:1234
## Median : 849.0 Median :1566 Median :1770 Median :1766 Median :1918
## Mean : 979.8 Mean :1561 Mean :1737 Mean :1732 Mean :1869
## 3rd Qu.:1241.5 3rd Qu.:2002 3rd Qu.:2242 3rd Qu.:2275 3rd Qu.:2448
## Max. :6146.0 Max. :6331 Max. :6394 Max. :6354 Max. :6308
##
## sentinel_o_b09 sentinel_o_b10
## Min. : 53.0 Min. : 28.0
## 1st Qu.: 600.5 1st Qu.: 309.2
## Median :1121.0 Median : 639.0
## Mean :1300.0 Mean : 834.0
## 3rd Qu.:1829.5 3rd Qu.:1186.0
## Max. :5973.0 Max. :4484.0
## En este caso se observa como no hay datos ausentes (NA). En caso de que aparecieran es importante eliminarlos, empleando para ello la función na.omit(). Una vez aplicada podemos emplear la función complete.cases() para comprobar que se han borrado.
train_data@data= na.omit(train_data@data)
complete.cases(train_data@data)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [113] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [127] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [141] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [155] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [183] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [197] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [225] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [239] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [253] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [267] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [281] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [295] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [309] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [323] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [337] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [351] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [365] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [379] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [393] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [407] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [435] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [449] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [463] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [477] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [491] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [505] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [519] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [533] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [547] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [561] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [575] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [589] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [603] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [617] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [631] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [645] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [659] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [673] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [687] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [701] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [715] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [729] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [743] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [757] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [771] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [785] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [799] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [813] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [827] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [841] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [855] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [869] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [883] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [897] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [911] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [925] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [939] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [953] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [967] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [981] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [995] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1009] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1023] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1037] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1051] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1065] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1079] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1093] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEPreparación del set de entrenamiento y testeo
Es recomendable separar de forma aleatoria el set de entrenamiento inicial en 3 grupos: entrenamiento, validación y testeo. En este caso solo lo vamos a separar en los dos primeros. Para ello, en primer lugar es necesario establecer un valor predefinido de semilla empleando set.seed().
hre_seed<- 123
set.seed(hre_seed)Ahora, dividiremos nuestro set de entrenamiento inicial en entrenamiento y testeo usando para ello la función createDataPartition() del paquete caret. Por ejemplo, vamos a establecer que el 80% de los datos iniciales pasen a ser de entrenamiento y el 20% de test. Recordar que el entrenamiento nos permitirá optimizar los parámetros iniciales del modelo mientras que los de testeo nos permitirán evaluar la calidad del mismo.
Nota: Se ha establecido como variable train_data@data$leyenda pues es esta la que que contiene la etiqueta de las clases. Por otro lado el parámetro p contiene el porcentaje a emplear en la separación de la muestra. Finalmente, el parámetro list indica si devuelve una lista o una matriz, en nuestro caso indicaremos FALSE de forma que devuelve una matriz.
Así, el resultado de la variable inTraining como podemos comprobar es una lista de valores núméricos indicando el índice de los elementos empleados para entrenamiento.
inTraining <- createDataPartition(train_data@data$leyenda, p=0.80,list=FALSE)
training <-train_data@data[inTraining,]
training=training[,-(1:4)] #Borramos columnas no necesarias
testing <- train_data@data[-inTraining,]
testing=testing[,-(1:4)] #Borramos columnas no necesariasResumen estadistico de los set de entreneamiento y testeo
Antes de comenzar el entrenamiento del clasificador de machine learning previo a la clasificación es necesario realizar un chequeo de los set de datos pues puede ser que las imágenes presenten problemas o que hayamos cometido errores en la identificación. Así, en primer lugar vamos a obtener un resumen estadístico de ambos set.
summary(training)## leyenda sentinel_o_b01 sentinel_o_b02 sentinel_o_b03
## Pinos :160 Min. : 761.0 Min. : 457.0 Min. : 249.0
## Bosques Mixtos:120 1st Qu.: 884.0 1st Qu.: 624.0 1st Qu.: 402.0
## Veg. Dispersa :120 Median : 983.5 Median : 792.0 Median : 593.5
## Matorral : 80 Mean :1100.7 Mean : 903.2 Mean : 775.1
## Quercineas : 80 3rd Qu.:1126.0 3rd Qu.: 995.8 3rd Qu.: 973.8
## -- : 40 Max. :5342.0 Max. :5290.0 Max. :5925.0
## (Other) :280
## sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## Min. : 226.0 Min. : 247 Min. : 233 Min. : 212 Min. : 187
## 1st Qu.: 543.0 1st Qu.:1017 1st Qu.:1137 1st Qu.:1114 1st Qu.:1204
## Median : 862.0 Median :1560 Median :1771 Median :1774 Median :1924
## Mean : 978.8 Mean :1555 Mean :1729 Mean :1723 Mean :1862
## 3rd Qu.:1243.0 3rd Qu.:2010 3rd Qu.:2250 3rd Qu.:2278 3rd Qu.:2443
## Max. :6146.0 Max. :6331 Max. :6394 Max. :6354 Max. :6308
##
## sentinel_o_b09 sentinel_o_b10
## Min. : 53 Min. : 28.0
## 1st Qu.: 580 1st Qu.: 301.0
## Median :1138 Median : 646.5
## Mean :1315 Mean : 845.8
## 3rd Qu.:1848 3rd Qu.:1196.5
## Max. :5973 Max. :4484.0
## summary(testing)## leyenda sentinel_o_b01 sentinel_o_b02 sentinel_o_b03
## Pinos :40 Min. : 783.0 Min. : 496.0 Min. : 282.0
## Bosques Mixtos:30 1st Qu.: 899.5 1st Qu.: 662.0 1st Qu.: 434.0
## Veg. Dispersa :30 Median : 972.0 Median : 772.5 Median : 568.0
## Matorral :20 Mean :1110.8 Mean : 912.7 Mean : 775.9
## Quercineas :20 3rd Qu.:1123.8 3rd Qu.: 995.2 3rd Qu.: 917.5
## -- :10 Max. :5384.0 Max. :5289.0 Max. :5881.0
## (Other) :70
## sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## Min. : 264.0 Min. : 303 Min. : 267 Min. : 253 Min. : 213
## 1st Qu.: 612.5 1st Qu.:1154 1st Qu.:1264 1st Qu.:1254 1st Qu.:1335
## Median : 807.0 Median :1576 Median :1751 Median :1749 Median :1904
## Mean : 983.7 Mean :1588 Mean :1766 Mean :1767 Mean :1900
## 3rd Qu.:1219.2 3rd Qu.:1958 3rd Qu.:2195 3rd Qu.:2262 3rd Qu.:2457
## Max. :5873.0 Max. :5905 Max. :5943 Max. :6284 Max. :6112
##
## sentinel_o_b09 sentinel_o_b10
## Min. : 87.0 Min. : 47.0
## 1st Qu.: 645.8 1st Qu.: 320.0
## Median :1094.0 Median : 596.5
## Mean :1241.2 Mean : 787.0
## 3rd Qu.:1681.2 3rd Qu.:1072.5
## Max. :4716.0 Max. :3370.0
## Posteriormente vamos a generar un grafico de densidades para cada clase / banda que permita representar la distribución de los datos. Esto va a permitirnos evaluar si hay una adecuada separabilidad entre las clases. Además permite determinar si el efecto cizalla en la distribución es acusado o no. Nota: Se han seleccionado los índices 2 al 11 pues en el caso de este ejemplo contienen los datos de reflectancia para cada punto en todas las bandas empleadas.
featurePlot(x=training[,2:11],
y=training$leyenda,
plot="density",
labels=c("Reflectancia","Distribucion densidades"),
layout=c(2,2))
Por otro lado podemos calcular la cizalla mediante la función skewness() de la librería e1071. Nos apartorá información si al distribución es simétrica o no. Por lo general, una distribución es simétrica cuando el valor de skewness es 0 o próximo a 0.
skwenessvalues <- apply(training[,2:11],2,skewness)
skwenessvalues## sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05
## 4.6181023 3.6657509 3.0889411 2.4448703 0.8939717
## sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10
## 0.5272516 0.4532052 0.3274439 1.1063450 1.4029482Por otro lado, si se detecta alguna distribución bimodal puede ser indicativo de una posible presencia de errores groseros en el muestreo. De forma complementaria podemos representar los datos mediante cajas de bigotes con objeto de ver esta presencia.
featurePlot(x=training[,2:11],
y=training$leyenda,
plot="box",
layout=c(2,2))
La posible presencia de errores groseros podría deberse por una parte a fallos humanos o tambien a la propia variabilidad del territorio y el comportamiento de las clases. Supongamos por ejemplo que contamos con las clases “uso agricola” y “suelo desnudo”, es posible que estas dos clases presenten un comportamiento similar, siendo adecuado analizar la correlación entre bandas.
A modo de ejemplo se presentan graficos de correlación entre dos clases y 6 bandas espectrales, viendo una clara correlación entre bandas.
band1_2 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b02))+
geom_point(aes(shape=leyenda,colour=leyenda))
band1_3 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b03))+
geom_point(aes(shape=leyenda,colour=leyenda))
band1_4 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b04))+
geom_point(aes(shape=leyenda,colour=leyenda))
band1_5 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b05))+
geom_point(aes(shape=leyenda,colour=leyenda))
grid.arrange(band1_2,band1_3,band1_4,band1_5)## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.## Warning: Removed 520 rows containing missing values (geom_point).## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.## Warning: Removed 520 rows containing missing values (geom_point).## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.## Warning: Removed 520 rows containing missing values (geom_point).## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.## Warning: Removed 520 rows containing missing values (geom_point).
Numéricamente, mediante la función cor() calculamos la correlacion entre las bandas espectrales de la escena. El resultado puede ser “complejo” de analizar, siendo mejor una representación gráfica.
bandcorrelaciones = cor(training[,2:11])
bandcorrelaciones## sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04
## sentinel_o_b01 1.0000000 0.9859977 0.9640630 0.9285340
## sentinel_o_b02 0.9859977 1.0000000 0.9882650 0.9734298
## sentinel_o_b03 0.9640630 0.9882650 1.0000000 0.9811764
## sentinel_o_b04 0.9285340 0.9734298 0.9811764 1.0000000
## sentinel_o_b05 0.7511718 0.8281299 0.8091457 0.8948271
## sentinel_o_b06 0.6895358 0.7743948 0.7554891 0.8514766
## sentinel_o_b07 0.6781686 0.7679473 0.7539346 0.8491108
## sentinel_o_b08 0.6479488 0.7417461 0.7300217 0.8323305
## sentinel_o_b09 0.4105274 0.5346405 0.6079181 0.6789671
## sentinel_o_b10 0.4526035 0.5684750 0.6488496 0.6984570
## sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## sentinel_o_b01 0.7511718 0.6895358 0.6781686 0.6479488
## sentinel_o_b02 0.8281299 0.7743948 0.7679473 0.7417461
## sentinel_o_b03 0.8091457 0.7554891 0.7539346 0.7300217
## sentinel_o_b04 0.8948271 0.8514766 0.8491108 0.8323305
## sentinel_o_b05 1.0000000 0.9912350 0.9877215 0.9837551
## sentinel_o_b06 0.9912350 1.0000000 0.9936050 0.9932052
## sentinel_o_b07 0.9877215 0.9936050 1.0000000 0.9938419
## sentinel_o_b08 0.9837551 0.9932052 0.9938419 1.0000000
## sentinel_o_b09 0.6890304 0.6935422 0.7193895 0.7380472
## sentinel_o_b10 0.6491981 0.6426895 0.6703054 0.6839346
## sentinel_o_b09 sentinel_o_b10
## sentinel_o_b01 0.4105274 0.4526035
## sentinel_o_b02 0.5346405 0.5684750
## sentinel_o_b03 0.6079181 0.6488496
## sentinel_o_b04 0.6789671 0.6984570
## sentinel_o_b05 0.6890304 0.6491981
## sentinel_o_b06 0.6935422 0.6426895
## sentinel_o_b07 0.7193895 0.6703054
## sentinel_o_b08 0.7380472 0.6839346
## sentinel_o_b09 1.0000000 0.9862417
## sentinel_o_b10 0.9862417 1.0000000Esta sería la representación gráfica de la matriz de correlación
corrplot(bandcorrelaciones,method="number")
corrplot(bandcorrelaciones,method="number",type = "upper")
corrplot(bandcorrelaciones,method="color",type="lower")
Definición de los parámetros del modelo para entrenamiento
Este paso es uno de los mas importantes, pued de la correcta configuración de los parámetros dependerán nuestros resultados. Para ello, usaremos la función trainControl() dentro del paquete caret. Esta se va a encargar de definir la configuración óptima del modelo. La función presenta tres parámetros:
method: “boot”, “boot632”, “optimism_boot”, “boot_all”, “cv”, “repeatedcv”, “LOOCV”,etc…
number: establece el numero de partes o bloques a dividir el conjunto de datos del mismo tamaño.
repeat: número de repeticiones.
La función selecciona el valor que da el mejor resultado.
fitControl <- trainControl(method = "repeatedcv",
number=5,
repeats=5)Entrenamiento de un clasificador KNN (k-nearest neighbors)
La función train() del paquete caret es la que contiene la lógica para definir el clasificador de machine learning. Aspectos a destacar en el uso de la función:
clases ~ : Indica que se emplearán todos los atributos.
El parámetro data contiene las variables predictoras.
El parámetro method indica el métdo del clasificador a emplear, en este caso knn.
trControl establece como vamos a definir los parámetros del modelo, ver paso anterior.
En el caso de usar un clasificador KNN como método no paramétrico utilizado en clasificación, este predice o clasifica una muestra de datos utilizando la muestra más cercana a k de los datos de entrenamiento.
Por tanto, un clasificador KNN depende en gran medida de cómo se defina la distancia entre las muestras. Aunque hay muchas métricas de distancia, la distancia euclidiana es la que se utiliza habitualmente.
Dado que la distancia entre las muestras es crítica, se recomienda preprocesar (centrar y escalar) las variables predictoras antes de ejecutar el clasificador KNN. Esto elimina los sesgos y permite que todos los predictores sean tratados por igual al calcular la distancia.
Las ventajas del clasificador KNN son:
Es un clasificador sencillo que puede implementarse fácilmente.
Resulta apropiado para manejar clases multimodales.
Como inconveniente, requiere el cálculo de la distancia de los vecinos más cercanos, lo que puede demandar una alta capacidad de cómputo, sobre todo si el conjunto de datos de entrenamiento es grande.
knnFit <- train(leyenda ~ .,data=training,
method="kknn",
preProcess=c("center","scale"),
trControl = fitControl)print (knnFit)## k-Nearest Neighbors
##
## 880 samples
## 10 predictor
## 13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa'
##
## Pre-processing: centered (10), scaled (10)
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ...
## Resampling results across tuning parameters:
##
## kmax Accuracy Kappa
## 5 0.3715909 0.2954465
## 7 0.3804545 0.3045002
## 9 0.3822727 0.3054572
##
## Tuning parameter 'distance' was held constant at a value of 2
## Tuning
## parameter 'kernel' was held constant at a value of optimal
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were kmax = 9, distance = 2 and kernel
## = optimal.Si graficamos los resultados podemos observar cual es el número adecuado de vecinos.
plot (knnFit)
Así como la configuaración del modelo
knnFit$finalModel##
## Call:
## kknn::train.kknn(formula = .outcome ~ ., data = dat, kmax = param$kmax, distance = param$distance, kernel = as.character(param$kernel))
##
## Type of response variable: nominal
## Minimal misclassification: 0.6079545
## Best kernel: optimal
## Best k: 7Y la importancia de las variables empleadas a través de la función varImp().
knnvarImp <-varImp(knnFit,compete=FALSE)
plot(knnvarImp)
Además del entrenamiento, resulta necesario realizar un testeo del mismo, empleando para ello la función predict()
pred_knnFit <-predict(knnFit,newdata=testing)Una vez realizado la fase de testeo sera necesario realizar un análisis de la calidad obtenida
confusionMatrix(data=pred_knnFit,testing$leyenda)## Confusion Matrix and Statistics
##
## Reference
## Prediction -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 3 2 0 0
## Bosques Mixtos 1 14 1 0
## Cadufolios 0 2 2 1
## Castaños/Caducifolios 0 0 0 9
## Eucalipto 1 2 0 0
## Matorral 1 3 1 0
## Pinos 4 3 0 0
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 1 0 0
## Quercineas 0 0 0 0
## Ribera 0 1 0 0
## Suelos? 0 2 1 0
## Veg. Dispersa 0 0 5 0
## Reference
## Prediction Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 0 1 1 0 0
## Bosques Mixtos 1 2 4 1 0
## Cadufolios 0 2 1 0 1
## Castaños/Caducifolios 0 0 0 0 0
## Eucalipto 4 0 1 1 0
## Matorral 1 6 1 0 0
## Pinos 2 3 20 4 3
## Pinos/Pinsapos 1 0 1 0 0
## Pinsapos 0 0 3 0 5
## Quercineas 0 0 1 4 1
## Ribera 1 0 2 0 0
## Suelos? 0 2 1 0 0
## Veg. Dispersa 0 4 4 0 0
## Reference
## Prediction Quercineas Ribera Suelos? Veg. Dispersa
## -- 0 1 1 0
## Bosques Mixtos 4 1 0 2
## Cadufolios 0 0 0 2
## Castaños/Caducifolios 0 0 0 0
## Eucalipto 0 0 1 1
## Matorral 2 0 2 7
## Pinos 2 2 2 1
## Pinos/Pinsapos 1 1 0 1
## Pinsapos 0 1 0 0
## Quercineas 9 2 1 0
## Ribera 0 1 0 0
## Suelos? 0 0 0 2
## Veg. Dispersa 2 1 3 14
##
## Overall Statistics
##
## Accuracy : 0.3955
## 95% CI : (0.3304, 0.4634)
## No Information Rate : 0.1818
## P-Value [Acc > NIR] : 1.086e-13
##
## Kappa : 0.3214
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity 0.30000 0.46667 0.200000
## Specificity 0.97143 0.91053 0.957143
## Pos Pred Value 0.33333 0.45161 0.181818
## Neg Pred Value 0.96682 0.91534 0.961722
## Prevalence 0.04545 0.13636 0.045455
## Detection Rate 0.01364 0.06364 0.009091
## Detection Prevalence 0.04091 0.14091 0.050000
## Balanced Accuracy 0.63571 0.68860 0.578571
## Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity 0.90000 0.40000
## Specificity 1.00000 0.96667
## Pos Pred Value 1.00000 0.36364
## Neg Pred Value 0.99526 0.97129
## Prevalence 0.04545 0.04545
## Detection Rate 0.04091 0.01818
## Detection Prevalence 0.04091 0.05000
## Balanced Accuracy 0.95000 0.68333
## Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity 0.30000 0.50000 0.00000
## Specificity 0.91000 0.85556 0.97619
## Pos Pred Value 0.25000 0.43478 0.00000
## Neg Pred Value 0.92857 0.88506 0.95349
## Prevalence 0.09091 0.18182 0.04545
## Detection Rate 0.02727 0.09091 0.00000
## Detection Prevalence 0.10909 0.20909 0.02273
## Balanced Accuracy 0.60500 0.67778 0.48810
## Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity 0.50000 0.45000 0.100000
## Specificity 0.97619 0.95500 0.980952
## Pos Pred Value 0.50000 0.50000 0.200000
## Neg Pred Value 0.97619 0.94554 0.958140
## Prevalence 0.04545 0.09091 0.045455
## Detection Rate 0.02273 0.04091 0.004545
## Detection Prevalence 0.04545 0.08182 0.022727
## Balanced Accuracy 0.73810 0.70250 0.540476
## Class: Suelos? Class: Veg. Dispersa
## Sensitivity 0.00000 0.46667
## Specificity 0.96190 0.90000
## Pos Pred Value 0.00000 0.42424
## Neg Pred Value 0.95283 0.91444
## Prevalence 0.04545 0.13636
## Detection Rate 0.00000 0.06364
## Detection Prevalence 0.03636 0.15000
## Balanced Accuracy 0.48095 0.68333Si los resultados en el control de calidad son positivos quedaría aplicar el modelo sobre la imagen para obtener la clasificación. (Este proceso requiere tiempo de computo dependiendo de la máquina empleada)
LC_knnFit <- predict(sentinel_o,knnFit)mapview(LC_knnFit)## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 4220964
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 4220964 '## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definitionEntrenar a un clasificador ANN (Artificial Neural Networks)
Un clasificador ANN es un método basado en simular el funcionamiento del cerebro humano para: a) adquisición de conocimientos, b) recordar, c) sintetizar y, d) resolver problemas. Hay distintos clasificadore KNN com el perceptrón multicapa (MLP); mapas de características auto-organizados (SOM) de Kohonen; las redes de Hopfield; el clasificador de Carpenter/Grossberg.
De todos ellos, el de tipo MLP es uno de los modelos de red neuronal más utilizados, generalmente consta de tres o más capas:
Una capa de entrada: consta de uno o varios elementos de procesamiento que presentan los datos de entrenamiento.
Una o más capas ocultas, responsable de la representación interna de los datos, así como de la transformación de la información entre las capas de entrada y de salida.
Una capa de salida: consta de uno o varios elementos de procesamiento que almacenan los resultados de la red.
Las ventajas de un clasificadores KNN son que:
Los datos no deben estar sujetos a seguir una distribución normal.
Tienen la capacidad de generar límites de decisión no lineales.
Capacidad de aprender patrones complejos.
Sin embargo, los clasificadores de KNN son propensos al overfiting sobreajuste, siendo complejo el diseñar una red neuronal eficaz. El overfitting se debe a que el modelo se ajustará a aprender los casos particulares que le enseñamos, siendo incapaz de reconocer nuevos datos de entrada.
annFit <- train(leyenda ~ ., data = training,
method = "nnet",
preProcess = c("center", "scale"),
trControl = fitControl)print (annFit)## Neural Network
##
## 880 samples
## 10 predictor
## 13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa'
##
## Pre-processing: centered (10), scaled (10)
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ...
## Resampling results across tuning parameters:
##
## size decay Accuracy Kappa
## 1 0e+00 0.2854545 0.1586734
## 1 1e-04 0.2815909 0.1543827
## 1 1e-01 0.2643182 0.1236267
## 3 0e+00 0.3820455 0.2946602
## 3 1e-04 0.3943182 0.3087801
## 3 1e-01 0.3850000 0.2908279
## 5 0e+00 0.4190909 0.3427078
## 5 1e-04 0.4245455 0.3486156
## 5 1e-01 0.4234091 0.3399757
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were size = 5 and decay = 1e-04.La siguiente figura muestra la relación entre pesos y capas ocultas de la red. En este ejemplo la mejor configuración del modelo ANN se obtiene con 5 nodos y un peso igual a 0.1.
plot (annFit)
Así, la configuración del modelo ANN es: - Se ha entrenado con 10 variables (las 10 bandas espectrales de Sentinel 2) - 5 nodos en la capa oculta de la red - 9 clase para la capa de salida
annFit$finalModel## a 10-5-13 network with 133 weights
## inputs: sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10
## output(s): .outcome
## options were - softmax modelling decay=1e-04Si deseamos visualizar la red emplearemos la función plotnet.
plotnet(annFit$finalModel)
Y mediante la función olden podemos analizar a importancia relativa de las variables predictoras. En este caso la banda 1 es la que más importancia presenta, todo lo contrario que la banda 5.
olden(annFit)
Y ahora, realizamos la predicción y el control de calidad sobre esta.
pred_annFit<- predict(annFit, newdata = testing)
confusionMatrix(data = pred_annFit, testing$leyenda)## Confusion Matrix and Statistics
##
## Reference
## Prediction -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 2 0 0 0
## Bosques Mixtos 2 15 0 0
## Cadufolios 0 1 1 0
## Castaños/Caducifolios 0 0 0 9
## Eucalipto 0 1 0 0
## Matorral 0 1 2 0
## Pinos 4 7 1 0
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 0 0 0
## Quercineas 0 1 1 0
## Ribera 0 0 1 0
## Suelos? 0 0 0 0
## Veg. Dispersa 2 4 4 1
## Reference
## Prediction Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 0 3 1 0 0
## Bosques Mixtos 4 6 6 3 0
## Cadufolios 0 0 1 0 1
## Castaños/Caducifolios 0 0 0 0 0
## Eucalipto 5 0 1 0 0
## Matorral 0 6 0 0 0
## Pinos 0 1 28 6 5
## Pinos/Pinsapos 0 0 0 0 0
## Pinsapos 0 2 0 0 2
## Quercineas 0 0 0 0 0
## Ribera 1 0 0 0 0
## Suelos? 0 0 0 0 0
## Veg. Dispersa 0 2 3 1 2
## Reference
## Prediction Quercineas Ribera Suelos? Veg. Dispersa
## -- 1 0 0 2
## Bosques Mixtos 3 4 1 2
## Cadufolios 0 0 0 0
## Castaños/Caducifolios 0 0 0 1
## Eucalipto 1 1 0 0
## Matorral 0 0 2 8
## Pinos 1 3 5 2
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 0 0 0
## Quercineas 10 1 0 3
## Ribera 0 0 0 0
## Suelos? 0 0 0 2
## Veg. Dispersa 4 1 2 10
##
## Overall Statistics
##
## Accuracy : 0.4
## 95% CI : (0.3347, 0.468)
## No Information Rate : 0.1818
## P-Value [Acc > NIR] : 3.613e-14
##
## Kappa : 0.3138
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity 0.200000 0.50000 0.100000
## Specificity 0.966667 0.83684 0.985714
## Pos Pred Value 0.222222 0.32609 0.250000
## Neg Pred Value 0.962085 0.91379 0.958333
## Prevalence 0.045455 0.13636 0.045455
## Detection Rate 0.009091 0.06818 0.004545
## Detection Prevalence 0.040909 0.20909 0.018182
## Balanced Accuracy 0.583333 0.66842 0.542857
## Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity 0.90000 0.50000
## Specificity 0.99524 0.98095
## Pos Pred Value 0.90000 0.55556
## Neg Pred Value 0.99524 0.97630
## Prevalence 0.04545 0.04545
## Detection Rate 0.04091 0.02273
## Detection Prevalence 0.04545 0.04091
## Balanced Accuracy 0.94762 0.74048
## Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity 0.30000 0.7000 0.00000
## Specificity 0.93500 0.8056 1.00000
## Pos Pred Value 0.31579 0.4444 NaN
## Neg Pred Value 0.93035 0.9236 0.95455
## Prevalence 0.09091 0.1818 0.04545
## Detection Rate 0.02727 0.1273 0.00000
## Detection Prevalence 0.08636 0.2864 0.00000
## Balanced Accuracy 0.61750 0.7528 0.50000
## Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity 0.200000 0.50000 0.000000
## Specificity 0.990476 0.97000 0.990476
## Pos Pred Value 0.500000 0.62500 0.000000
## Neg Pred Value 0.962963 0.95098 0.954128
## Prevalence 0.045455 0.09091 0.045455
## Detection Rate 0.009091 0.04545 0.000000
## Detection Prevalence 0.018182 0.07273 0.009091
## Balanced Accuracy 0.595238 0.73500 0.495238
## Class: Suelos? Class: Veg. Dispersa
## Sensitivity 0.000000 0.33333
## Specificity 0.990476 0.86316
## Pos Pred Value 0.000000 0.27778
## Neg Pred Value 0.954128 0.89130
## Prevalence 0.045455 0.13636
## Detection Rate 0.000000 0.04545
## Detection Prevalence 0.009091 0.16364
## Balanced Accuracy 0.495238 0.59825Entrenar a un clasificador SVM (Support Vector Machine)
Un clasificador SVM se basa en modelos estadísticos los cuales, con caracter tiene por objeto localizar una frontera de decisión (separación) óptima que maximice el margen entre dos clases.
La ubicación del límite de decisión dependerá solo de un subconjunto de puntos de datos de entrenamiento que están más cerca de él. Este subconjunto de puntos de datos de entrenamiento más cercanos al límite de decisión se conoce como vectores de soporte. A lo largo del tiempo han aparecido evolutivos incorporando funciones de coste, funciones para permitir límites de clase no lineales.
Además, funciones polinómicas, base radial o tangente hiperbólica, se desarrollaron para transformar el conjunto de datos de entrenamiento en un espacio de características de mayor dimensión para los problemas de clasificación no lineal.
Las ventajas de un clasificador SVM son:
Se elimina el problema de la posible suposición de distribución normal de los datos.
Uso de una funciones para resolver problemas complejos.
Se adapta relativamente bien a datos de alta dimensión.
No obstante, como inconventiente, este clasificador necesita de más tiempo de entrenamiento.
svm_model<-train(leyenda~.,data=training,
method = "svmRadial",
trControl = fitControl,
preProc = c("center", "scale"),
tuneLength = 3)Al hacer una impresión del modelo se observa como este tiene una calidad igual a 0.36 y un coste igual a 1.
print (svm_model)## Support Vector Machines with Radial Basis Function Kernel
##
## 880 samples
## 10 predictor
## 13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa'
##
## Pre-processing: centered (10), scaled (10)
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ...
## Resampling results across tuning parameters:
##
## C Accuracy Kappa
## 0.25 0.3131818 0.1843286
## 0.50 0.3295455 0.2112172
## 1.00 0.3665909 0.2671307
##
## Tuning parameter 'sigma' was held constant at a value of 0.5696137
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.5696137 and C = 1.plot (svm_model)
El siguiente paso es mostrar la importancia de las variables de entrada
svm_varImp <- varImp(svm_model, compete = FALSE)
plot(svm_varImp)
Y finalmente realizar un control de calidad
pred_svm<- predict(svm_model, newdata = testing)
confusionMatrix(data = pred_svm, testing$leyenda)## Confusion Matrix and Statistics
##
## Reference
## Prediction -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 3 0 0 0
## Bosques Mixtos 1 5 0 0
## Cadufolios 0 2 0 2
## Castaños/Caducifolios 0 0 0 8
## Eucalipto 0 0 1 0
## Matorral 0 3 2 0
## Pinos 6 16 1 0
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 1 1 0
## Quercineas 0 2 0 0
## Ribera 0 0 0 0
## Suelos? 0 0 0 0
## Veg. Dispersa 0 1 5 0
## Reference
## Prediction Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 0 2 0 0 0
## Bosques Mixtos 1 4 2 2 0
## Cadufolios 0 1 0 0 0
## Castaños/Caducifolios 0 0 0 0 0
## Eucalipto 0 0 1 0 0
## Matorral 0 10 0 0 1
## Pinos 6 2 33 8 6
## Pinos/Pinsapos 0 0 0 0 0
## Pinsapos 2 0 1 0 2
## Quercineas 0 0 0 0 1
## Ribera 0 0 0 0 0
## Suelos? 0 0 0 0 0
## Veg. Dispersa 1 1 3 0 0
## Reference
## Prediction Quercineas Ribera Suelos? Veg. Dispersa
## -- 0 0 0 0
## Bosques Mixtos 2 0 1 3
## Cadufolios 0 0 0 0
## Castaños/Caducifolios 0 0 0 0
## Eucalipto 0 0 1 0
## Matorral 0 0 1 5
## Pinos 8 7 3 2
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 0 0 0
## Quercineas 7 3 0 2
## Ribera 0 0 0 0
## Suelos? 0 0 0 0
## Veg. Dispersa 3 0 4 18
##
## Overall Statistics
##
## Accuracy : 0.3909
## 95% CI : (0.326, 0.4588)
## No Information Rate : 0.1818
## P-Value [Acc > NIR] : 3.202e-13
##
## Kappa : 0.2939
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity 0.30000 0.16667 0.00000
## Specificity 0.99048 0.91579 0.97619
## Pos Pred Value 0.60000 0.23810 0.00000
## Neg Pred Value 0.96744 0.87437 0.95349
## Prevalence 0.04545 0.13636 0.04545
## Detection Rate 0.01364 0.02273 0.00000
## Detection Prevalence 0.02273 0.09545 0.02273
## Balanced Accuracy 0.64524 0.54123 0.48810
## Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity 0.80000 0.00000
## Specificity 1.00000 0.98571
## Pos Pred Value 1.00000 0.00000
## Neg Pred Value 0.99057 0.95392
## Prevalence 0.04545 0.04545
## Detection Rate 0.03636 0.00000
## Detection Prevalence 0.03636 0.01364
## Balanced Accuracy 0.90000 0.49286
## Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity 0.50000 0.8250 0.00000
## Specificity 0.94000 0.6389 1.00000
## Pos Pred Value 0.45455 0.3367 NaN
## Neg Pred Value 0.94949 0.9426 0.95455
## Prevalence 0.09091 0.1818 0.04545
## Detection Rate 0.04545 0.1500 0.00000
## Detection Prevalence 0.10000 0.4455 0.00000
## Balanced Accuracy 0.72000 0.7319 0.50000
## Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity 0.200000 0.35000 0.00000
## Specificity 0.976190 0.96000 1.00000
## Pos Pred Value 0.285714 0.46667 NaN
## Neg Pred Value 0.962441 0.93659 0.95455
## Prevalence 0.045455 0.09091 0.04545
## Detection Rate 0.009091 0.03182 0.00000
## Detection Prevalence 0.031818 0.06818 0.00000
## Balanced Accuracy 0.588095 0.65500 0.50000
## Class: Suelos? Class: Veg. Dispersa
## Sensitivity 0.00000 0.60000
## Specificity 1.00000 0.90526
## Pos Pred Value NaN 0.50000
## Neg Pred Value 0.95455 0.93478
## Prevalence 0.04545 0.13636
## Detection Rate 0.00000 0.08182
## Detection Prevalence 0.00000 0.16364
## Balanced Accuracy 0.50000 0.75263Entrenar a un clasificador RF (Random Forest)
El clasificador RF es un método de aprendizaje automático de conjunto, que utiliza el muestreo bootstrap para construir muchos modelos de árboles de decisión individuales.
Usa un subconjunto aleatorio de variables predictoras (por ejemplo, las bandas Sentinel) para dividir los datos de observación en subconjuntos homogéneos, que se utilizan para construir cada modelo de árbol de decisión y una predicción. Luego, se promedian las predicciones del modelo de árbol de decisión individual para producir el etiquetado final.
El clasificador RF se ha utilizado con éxito para la clasificación de imágenes de teledetección porque presentan estas ventajas:
Permiten manejar grandes cantidades de datos.
Están libres de supuestos de distribución normal.
Son robustos a los valores atípicos y al ruido
Sin embargo, como inconventientes:
No es fácil interpretar los resultados del modelo RF.
Esta sesgado a favor de las variables predictoras con muchos niveles de categorías diferentes.
rf_model<-train(leyenda~.,data=training, method="rf",
trControl=fitControl,
prox=TRUE,
fitBest = FALSE,
returnData = TRUE)print(rf_model)## Random Forest
##
## 880 samples
## 10 predictor
## 13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.4195455 0.3434671
## 6 0.4284091 0.3542228
## 10 0.4304545 0.3566526
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 10.plot(rf_model)
rf_model$finalModel##
## Call:
## randomForest(x = x, y = y, mtry = min(param$mtry, ncol(x)), proximity = TRUE, fitBest = FALSE, returnData = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 10
##
## OOB estimate of error rate: 56.82%
## Confusion matrix:
## -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 24 4 0 0
## Bosques Mixtos 7 51 1 0
## Cadufolios 0 1 16 1
## Castaños/Caducifolios 0 1 1 36
## Eucalipto 0 6 0 0
## Matorral 3 7 4 1
## Pinos 1 25 3 0
## Pinos/Pinsapos 1 6 1 0
## Pinsapos 0 2 0 0
## Quercineas 2 10 1 0
## Ribera 0 5 0 0
## Suelos? 0 5 4 0
## Veg. Dispersa 0 9 12 4
## Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 0 1 7 3 0
## Bosques Mixtos 1 4 26 0 3
## Cadufolios 2 2 3 0 0
## Castaños/Caducifolios 0 0 0 0 0
## Eucalipto 11 1 12 0 0
## Matorral 0 26 9 1 0
## Pinos 3 4 93 6 6
## Pinos/Pinsapos 0 1 15 5 4
## Pinsapos 1 1 11 1 20
## Quercineas 1 0 18 0 1
## Ribera 2 3 6 0 0
## Suelos? 3 6 4 1 3
## Veg. Dispersa 5 18 14 1 0
## Quercineas Ribera Suelos? Veg. Dispersa class.error
## -- 0 0 0 1 0.4000000
## Bosques Mixtos 12 0 2 13 0.5750000
## Cadufolios 0 0 0 15 0.6000000
## Castaños/Caducifolios 0 0 0 2 0.1000000
## Eucalipto 0 4 2 4 0.7250000
## Matorral 2 1 2 24 0.6750000
## Pinos 9 2 1 7 0.4187500
## Pinos/Pinsapos 3 0 0 4 0.8750000
## Pinsapos 0 1 2 1 0.5000000
## Quercineas 34 3 1 9 0.5750000
## Ribera 4 16 0 4 0.6000000
## Suelos? 0 2 1 11 0.9750000
## Veg. Dispersa 6 1 3 47 0.6083333rf_varImp <- varImp(rf_model, compete = FALSE)
plot(rf_varImp)
Como en los clasificadores estudiados anteriormente realizaremos un control de calidad.
pred_rf <- predict(rf_model$finalModel,
newdata = testing)confusionMatrix(data = pred_rf, testing$leyenda)## Confusion Matrix and Statistics
##
## Reference
## Prediction -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## -- 3 2 0 0
## Bosques Mixtos 1 13 0 0
## Cadufolios 0 2 3 1
## Castaños/Caducifolios 0 0 0 9
## Eucalipto 0 1 1 0
## Matorral 0 2 1 0
## Pinos 6 6 0 0
## Pinos/Pinsapos 0 0 0 0
## Pinsapos 0 1 0 0
## Quercineas 0 2 0 0
## Ribera 0 0 1 0
## Suelos? 0 0 0 0
## Veg. Dispersa 0 1 4 0
## Reference
## Prediction Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## -- 2 2 1 0 0
## Bosques Mixtos 1 2 1 1 0
## Cadufolios 0 2 1 0 0
## Castaños/Caducifolios 0 0 0 0 0
## Eucalipto 1 0 2 1 1
## Matorral 0 8 0 0 0
## Pinos 3 1 22 6 5
## Pinos/Pinsapos 0 1 1 0 0
## Pinsapos 0 0 2 1 2
## Quercineas 0 0 3 0 0
## Ribera 1 0 0 0 0
## Suelos? 0 0 0 0 1
## Veg. Dispersa 2 4 7 1 1
## Reference
## Prediction Quercineas Ribera Suelos? Veg. Dispersa
## -- 0 0 1 0
## Bosques Mixtos 6 3 0 1
## Cadufolios 0 0 1 2
## Castaños/Caducifolios 0 0 0 0
## Eucalipto 0 0 0 0
## Matorral 1 0 2 5
## Pinos 4 4 3 2
## Pinos/Pinsapos 0 0 0 1
## Pinsapos 0 0 0 0
## Quercineas 7 1 0 0
## Ribera 0 2 0 0
## Suelos? 0 0 1 3
## Veg. Dispersa 2 0 2 16
##
## Overall Statistics
##
## Accuracy : 0.3955
## 95% CI : (0.3304, 0.4634)
## No Information Rate : 0.1818
## P-Value [Acc > NIR] : 1.086e-13
##
## Kappa : 0.3138
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity 0.30000 0.43333 0.30000
## Specificity 0.96190 0.91579 0.95714
## Pos Pred Value 0.27273 0.44828 0.25000
## Neg Pred Value 0.96651 0.91099 0.96635
## Prevalence 0.04545 0.13636 0.04545
## Detection Rate 0.01364 0.05909 0.01364
## Detection Prevalence 0.05000 0.13182 0.05455
## Balanced Accuracy 0.63095 0.67456 0.62857
## Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity 0.90000 0.100000
## Specificity 1.00000 0.971429
## Pos Pred Value 1.00000 0.142857
## Neg Pred Value 0.99526 0.957746
## Prevalence 0.04545 0.045455
## Detection Rate 0.04091 0.004545
## Detection Prevalence 0.04091 0.031818
## Balanced Accuracy 0.95000 0.535714
## Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity 0.40000 0.5500 0.00000
## Specificity 0.94500 0.7778 0.98571
## Pos Pred Value 0.42105 0.3548 0.00000
## Neg Pred Value 0.94030 0.8861 0.95392
## Prevalence 0.09091 0.1818 0.04545
## Detection Rate 0.03636 0.1000 0.00000
## Detection Prevalence 0.08636 0.2818 0.01364
## Balanced Accuracy 0.67250 0.6639 0.49286
## Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity 0.200000 0.35000 0.200000
## Specificity 0.980952 0.97000 0.990476
## Pos Pred Value 0.333333 0.53846 0.500000
## Neg Pred Value 0.962617 0.93720 0.962963
## Prevalence 0.045455 0.09091 0.045455
## Detection Rate 0.009091 0.03182 0.009091
## Detection Prevalence 0.027273 0.05909 0.018182
## Balanced Accuracy 0.590476 0.66000 0.595238
## Class: Suelos? Class: Veg. Dispersa
## Sensitivity 0.100000 0.53333
## Specificity 0.980952 0.87368
## Pos Pred Value 0.200000 0.40000
## Neg Pred Value 0.958140 0.92222
## Prevalence 0.045455 0.13636
## Detection Rate 0.004545 0.07273
## Detection Prevalence 0.022727 0.18182
## Balanced Accuracy 0.540476 0.70351Comparacón de los clasificadores estudiados
La comparación se realizará siguiendo una correlación cruzada, empleando para esto la función resample().
resamps <- resamples(list(knn = knnFit,
ann = annFit,
svm = svm_model,
rf = rf_model))Representamos mediante un boxplot los valores de accuracy y kappa para realizar una evaluación de la calidad entre clasificadores empleados.
bwplot(resamps, layout = c(3, 1))
A continuación se realizará la predicción de cada uno de los clasificadores anteriormente desarrollados.
LC_knnFit <-predict(sentinel_o,knnFit)
LC_ann <-predict(sentinel_o,annFit)
LC_svm <-predict(sentinel_o,svm_model)
LC_rf <-predict(sentinel_o,rf_model)Clasificación digital empleando datos multitemporales
El objetivo es obtener una clasificación de la misma zona de estudio que en el tutorial 1 pero en este caso considerando datos multitemporales correspondientes a tres escenas de Sentinel 2: primeravera, verano y otoño.
Vamos en primer lugar a generar un rasterstack según fecha.
dir_in='./Material_practicas/Sentinel/'
archivos <-list.files(paste(dir_in,"o",sep="")
,pattern = ".tif",
full.names = TRUE)
sentinel_o <- stack(archivos)
archivos <-list.files(paste(dir_in,"p",sep="")
,pattern = ".tif",
full.names = TRUE)
sentinel_p <- stack(archivos)
archivos <-list.files(paste(dir_in,"v",sep="")
,pattern = ".tif",
full.names = TRUE)
sentinel_v <- stack(archivos)plotRGB(sentinel_o,r=7,g=2,b=3,stretch="lin")
plotRGB(sentinel_p,r=7,g=2,b=3,stretch="lin")
plotRGB(sentinel_v,r=7,g=2,b=3,stretch="lin")
Crearemos a continuación un rasterstack resultante de las escenas Sentinel 2 de las tres fechas, teniendo un total de 30 bandas espectrales (10 por cada una de las fechas)
sentinel <- stack(sentinel_o,sentinel_p,sentinel_v)A partir de aquí, el desarrollo del clasificador será exactamente igual que en el caso de trabajar con una imagen, salvo que el tiempo de procesado será mayor